home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / src / exampleCode / audio / reverb / reverb.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-02  |  56.9 KB  |  1,836 lines
  1. /*
  2.  * Copyright (C) 1994, Silicon Graphics, Inc.
  3.  * All Rights Reserved.
  4.  *
  5.  * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
  6.  * the contents of this file may not be disclosed to third parties, copied or
  7.  * duplicated in any form, in whole or in part, without the prior written
  8.  * permission of Silicon Graphics, Inc.
  9.  *
  10.  * RESTRICTED RIGHTS LEGEND:
  11.  * Use, duplication or disclosure by the Government is subject to restrictions
  12.  * as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
  13.  * and Computer Software clause at DFARS 252.227-7013, and/or in similar or
  14.  * successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
  15.  * rights reserved under the Copyright Laws of the United States.
  16.  */
  17.  
  18.  
  19. /*********************************************************************** 
  20.  *
  21.  *  Simulated reverberation for real-time audio input and output.  
  22.  *
  23.  *  The algorithms generate room ambience rather than a realistic room 
  24.  *  response.  For proper operation, ensure the input and output sampling 
  25.  *  rates are equal and neither changes over time.
  26.  *
  27.  *      The most advanced algorithm in this code is moorer2er19, which
  28.  *      consumes roughly 50% of a 100 MHz R4000.  
  29.  *
  30.  *      For the pronounced effect, try:
  31.  *          > reverb -type moorer2er19 -reverbTime 10
  32.  *
  33.  *      For reverberation times > 30 seconds, the decay seems to fizzle out 
  34.  *      quite abruptly.  The cause is likely caused by single-precision 
  35.  *      floating point artifacts in the computation of the low-pass filter 
  36.  *      nested in each comb filter.  
  37.  *
  38.  *      -noMonoProcess      specifies that a reverberator is computed for 
  39.  *                          each input channel.  Otherwise, all inputs are 
  40.  *                          summed into a single reverberator.  Algorithm 
  41.  *                          types 'comb' and 'allpass'  are processed 
  42.  *                          with a feedback gain of opposite polarity in
  43.  *                          each of two output channels. 
  44.  *
  45.  *      -reverbTime         default is 8 seconds.  Applies only to comb, 
  46.  *                          allpass, schroeder2, moorer, moorer2, moorer2er7,
  47.  *                          moorer2er19 algorithms.  Others have decay time 
  48.  *                          as specified in design.
  49.  *
  50.  *      -noPrime            prevents reciculated delay buffer lengths from 
  51.  *                          quantization to prime value.  The "priming" 
  52.  *                          seems to reduce ringing due to transients such 
  53.  *                          as snare drum strike and low-level whining in 
  54.  *                          the presence of input noise, so this parameter 
  55.  *                          is provided for comparison.
  56.  *
  57.  *      -type:              defaults is moorer2
  58.  *
  59.  *      comb, allpass       comb and all-pass Interpolated Infinite Impulse 
  60.  *                          Response (IIIR) filters.  
  61.  *      schroeder1          quint series all-pass filter
  62.  *      schroeder2          quad parallel comb into dual series all-pass 
  63.  *                          filters
  64.  *      chamberlin          quint series all-pass filter
  65.  *      moorer              Moorer's hex parallel comb filters into single 
  66.  *                          all-pass filter
  67.  *      moorer2             hex parallel filtered-comb into single all-pass
  68.  *                          filters
  69.  *      er7, er19           early reflections patterns w/7 and 19 taps, 
  70.  *                          respectively 
  71.  *      moorer2er7, moorer2er19         early reflections + hex parallel 
  72.  *                          filtered-comb into single all-pass filters.  
  73.  *                          Early reflections fed into comb network.  
  74.  *                          Scatter in 19 tap ER pattern yields algorithm 
  75.  *                          with significantly less ring excited by 
  76.  *                          transients such as snare drum strike.
  77.  *
  78.  *      To my ears, moorer2 rings less to transients than moorer2er7, 
  79.  *      moorer2er19.  However, moorer2er7 and moorer2er19 have smoother 
  80.  *      decay than moorer2.  Not sure if ringing is caused by an error 
  81.  *      in my implementation or by a comb-filter quality imposed by the 
  82.  *      early reflection pattern.   Consider a diffusor other than the 
  83.  *      early reflection pattern. 
  84.  *
  85.  *  The following references contain detailed discussions of the 
  86.  *  non-proprietary implementations found in the code below:
  87.  *
  88.  *  + Schroeder, M.R. and Logan, B.F., "'Colorless' Artificial Reverberation."
  89.  *    Journal of the Audio Engineering Society v. 9, n. 3: pp. 192-197, 1962
  90.  *
  91.  *  + Schroeder, M.R. , "Natural Sounding Artificial Reverberation." Journal 
  92.  *    of the Audio Engineering Society v. 10, n. 3: pp. 219-223, 1962
  93.  *
  94.  *  + Moorer, J.A., "About This Reverberation Business", Foundations of 
  95.  *    Computer Music, Roads and Strawn, ed.: pp. 605-639, 1985
  96.  *
  97.  *  + Chamberlin, Hal, Musical Applications of Microprocessors, 2nd Ed., 
  98.  *    Hayden Books, A Division of Howard Sams & Co., Indianapolis, Indiana, 
  99.  *    pp. 508-512, 1985
  100.  *
  101.  *
  102.  *                Written by Gints Klimanis 
  103.  *            Silicon Graphics Computer Systems
  104.  *                    1994
  105.  *********************************************************************** */
  106. #include <stdio.h>
  107. #include <audio.h>
  108. #include <math.h>
  109. #include <stdlib.h>
  110. #include <sigfpe.h>
  111. #include <sys/schedctl.h>
  112.  
  113. #define FALSE    0
  114. #define TRUE    1
  115. #define LEFT    0
  116. #define RIGHT    1
  117.  
  118. #define MAX_CHANNELS    2
  119.  
  120. #define BLOCK_SIZE        512
  121. #define MAX_DELAY_SIZE        11025
  122. #define MAX_DELAY_COUNT        30
  123.  
  124. #define SINGLE_ALL_PASS            0
  125. #define SINGLE_COMB            1
  126. #define QUINT_ALL_PASS            2
  127. #define QUINT_ALL_PASS_STEREO        3
  128. #define QUAD_COMB_DUAL_ALL_PASS        4
  129. #define HEX_COMB_ALL_PASS        5
  130.  
  131. #define SINGLE_ALL_PASS_LPF        10
  132. #define SINGLE_COMB_LPF            11
  133. #define HEX_COMB_ALL_PASS_LPF        12
  134. #define HEX_COMB_ALL_PASS_LPF_ER7   13
  135. #define HEX_COMB_ALL_PASS_LPF_ER19  14
  136. #define ER7                15
  137. #define ER19                16
  138.  
  139. double    samplingRate = 44100;
  140. float    effectLevel = 0.03;
  141. char    effectType = HEX_COMB_ALL_PASS_LPF;
  142. char    stereoize = TRUE;
  143. char    snapToPrime = TRUE;    /* seems to reduce some ringing in my tests */
  144. char    monoProcess = TRUE;
  145. float    reverberationTime = 8;
  146.  
  147. char    verboseOperation = FALSE;
  148.  
  149. typedef struct {
  150.     float    *delayLines[MAX_CHANNELS][MAX_DELAY_COUNT], 
  151.         *delayOutputs[MAX_CHANNELS][MAX_DELAY_COUNT];
  152.     float    inLineDelays[MAX_CHANNELS][MAX_DELAY_COUNT];
  153.     int        delayLinePtrs[MAX_CHANNELS][MAX_DELAY_COUNT], 
  154.         delayLineLengths[MAX_CHANNELS][MAX_DELAY_COUNT];
  155.     double    delayLineTimes[MAX_CHANNELS][MAX_DELAY_COUNT];
  156.     float    delayLoopGains[MAX_CHANNELS][MAX_DELAY_COUNT];
  157.     float    inLineDelayGains[MAX_CHANNELS][MAX_DELAY_COUNT];
  158.  
  159.     double    reverberationTime;
  160.  
  161.     float    effectLevel;
  162.  
  163.     char    stereoize;
  164.     char    snapToPrime;
  165. } Reverb;
  166.  
  167. /* high frequency damping coefficients for 25000 & 50000 samples/second */
  168. float    highFrequencyDampingCoeffs_25000Hz[6] = {
  169.     0.24, 0.26, 0.28, 0.29, 0.30, 0.32};
  170. float    highFrequencyDampingCoeffs_50000Hz[6] = {
  171.     0.46, 0.48, 0.50, 0.52, 0.53, 0.55};
  172.  
  173. /* 
  174.  * Manfred Schroeder's conditions for artificial reverberators 
  175.  * (lifted from paper):
  176.  *
  177.  * 1) The frequency response must be flat when measured with narrow bands 
  178.  *    of noise, with the bandwidths corresponding to that of the transients 
  179.  *    in the sound to be reverberated.  This condition is, of course, 
  180.  *    fulfilled by reverberators which have a flat response even for 
  181.  *    sinusoidal excitation.
  182.  *
  183.  * 2) The normal modes of the reverberator must overlap and cover the entire
  184.  *    audio frequency range.
  185.  *
  186.  * 3) The reverberation times of the individual modes must be equal or 
  187.  *    nearly equal so that different frequency components of the sound will 
  188.  *    decay with equal rates.
  189.  *
  190.  * 4) The echo density a short interval after shock excitation must be high
  191.  *    enough so that individual echos are not resolved by the ear.
  192.  *
  193.  * 5) The echo response must be free from periodicities (flutter echos)
  194.  *
  195.  * 6) The amplitude-frequency response must not exhibit any apparent 
  196.  *    periodicities.  Periodic or comblike frequency responses produce an 
  197.  *    unpleasant hollow, reedy or metallic sound quality and give the 
  198.  *    impression that the sound is transmitted through a hollow tube or 
  199.  *    barrel.
  200.  */
  201.  
  202. static char applicationUsage[] = 
  203. "\t-help\n\
  204. \t-verbose\n\
  205. \t-level        0..1\n\
  206. \t-noMonoProcess    each input processed\n\
  207. \t-reverbTime        seconds\n\
  208. \t-noPrime\n\
  209. \t-type            {comb, allpass,\n\
  210. \t            schroeder1, schroeder2,\n\
  211. \t            chamberlin,\n\
  212. \t            moorer, moorer2,\n\
  213. \t            er7, er19\n\
  214. \t            moorer2er7, moorer2er19}\n";
  215.  
  216.     long 
  217. GetHardwareInputRate();
  218.     static void
  219. ParseCommandLine(int argc, char **argv, char fileName[][200]);
  220.     static void
  221. ProcessParameters(char type, double samplingRate, int channels, Reverb *reverb);
  222.  
  223.     void
  224. IIIRFloatsLPF(float *input, float *output, int length, 
  225.             float *delayLine, int delayLength, int *delayLineIndex, float loopGain, 
  226.             float *z, float loopGain2, 
  227.                 int topology);
  228.     static void
  229. ConvertDelayLengths(int *lengths, double *times, double samplingRate,
  230.             int count, char snapToPrime);
  231.  
  232.     static void
  233. MultiTapDelayFloats(float *input, float *output, int length, 
  234.             float *delayLine, int delayLength,
  235.                 float *tapGains, int *tapIndices, int tapCount);
  236.     int 
  237. IsPrime(int number);
  238.     int 
  239. NearestPrime(int number, int lowerBound, int upperBound);
  240.     int 
  241. UpperPrime(int number, int upperBound);
  242.     int 
  243. LowerPrime(int number, int lowerBound);
  244.     int
  245. SecondsToSamples(double time, double samplingRate);
  246.     void 
  247. InterleaveNFloatsToShorts(float *inputs[], short *output, int inputLength, 
  248.                 int interleaveFactor, char saturate);
  249.     void
  250. DelayFloats(float *in, float *out, long length, 
  251.         float *delayBuffer, int delayLength, int *delayBufferPtr,
  252.             char useCircularBuffer);
  253.     void 
  254. CopyFloats(float *in, float *out, int length);
  255.     void 
  256. SumFloats(float *in1,  float *in2, float *out, int length);
  257.     void 
  258. ScaleFloats(float *in, float *out, int length, float scaleFactor);
  259.     void
  260. IIIRFloats(float *in, float *out, int length, 
  261.             float *delayLine, int delayLength, int *delayLineIndex, 
  262.             float loopGain, int topology);
  263.     void 
  264. SumNFloats(float *inputs[], float *output, int length, int inputChannels);
  265.     void 
  266. DeinterleaveShortsToNFloats(short *input, float *outputs[], int outputLength,
  267.                 int interleaveFactor);
  268.     void 
  269. SetFloats(float *in, int length, float value);
  270.  
  271. /*********************************************************************** 
  272.  * main()
  273.  *********************************************************************** */
  274. main(int argc, char **argv)
  275. {
  276. int        i, j;
  277. int        someInteger;
  278. char        fileName[3][200];        
  279.  
  280. ALport        audioInPort, audioOutPort;  
  281. ALconfig    aconfig;
  282.  
  283. short        shortBuffer[BLOCK_SIZE*MAX_CHANNELS];
  284. float        *systemInputs[MAX_CHANNELS], *outputs[MAX_CHANNELS];
  285. float        **in, **out;
  286.  
  287. Reverb        *reverb;
  288. float        *gains, *inLineGains;
  289. int        *lengths, *delayLinePtrs;
  290. int        tapCount;
  291.  
  292. int        channelCount, channel;
  293. int        channelsToProcess;
  294. int        delayIndex, delayChannel;
  295. int        iiirType;
  296. float        *inPtr;
  297.  
  298. float        *tmpBuffer[10];
  299. FILE        *inHandle, *outHandle;
  300. char        string[100];
  301. int        priority;
  302.  
  303. ParseCommandLine(argc, argv, fileName);
  304.  
  305. /* Prevent OS from intervening w/floating point underflows.  
  306.     Set underflowing values to zero */
  307. sigfpe_[_UNDERFL].repls = _ZERO;
  308. handle_sigfpes(_ON, _EN_UNDERFL, NULL, _ABORT_ON_ERROR, NULL);
  309. /* try to schedule at higher priority to prevent drop outs */
  310. schedctl(NDPRI, 0, NDPNORMMAX);
  311.  
  312. /* initialize audio system */        
  313. aconfig = ALnewconfig();
  314. channelCount = MAX_CHANNELS;
  315. ALsetqueuesize(aconfig, BLOCK_SIZE*10);
  316. ALsetwidth(aconfig, AL_SAMPLE_16);
  317. ALsetchannels(aconfig, channelCount);
  318.  
  319. /* open audio ports w/ALconfiguration */
  320. audioInPort = ALopenport("xxx", "r", aconfig);
  321. audioOutPort = ALopenport("xxx", "w", aconfig);
  322. if ((audioInPort == NULL)||(audioOutPort == NULL))
  323.     {
  324.     fprintf(stderr, "main(): failed to open enough audio ports\n");
  325.     exit(-1);
  326.     }
  327. ALfreeconfig(aconfig);
  328. samplingRate = (double) GetHardwareInputRate();
  329.  
  330.  
  331. /* allocate memory */
  332. tmpBuffer[0] = (float *) malloc(2*BLOCK_SIZE*sizeof(float));
  333. tmpBuffer[1] = (float *) malloc(2*BLOCK_SIZE*sizeof(float));
  334. reverb = (Reverb *) malloc(sizeof(Reverb));
  335. for (i = 0; i < MAX_CHANNELS; i++)
  336.     {
  337.     systemInputs[i]  = (float *) malloc(BLOCK_SIZE*sizeof(float));
  338.     outputs[i] = (float *) malloc(BLOCK_SIZE*sizeof(float));
  339.  
  340.     for (j = 0; j < MAX_DELAY_COUNT; j++)
  341.     {
  342.     reverb->delayOutputs[i][j] = (float *) malloc(BLOCK_SIZE*sizeof(float));
  343.     reverb->delayLines[i][j] = (float *) malloc(MAX_DELAY_SIZE*sizeof(float));
  344.  
  345.     /* clear delay lines */
  346.     reverb->delayLinePtrs[i][j] = 0;
  347.     SetFloats(reverb->delayLines[i][j], MAX_DELAY_SIZE, 0);
  348.     reverb->inLineDelays[i][j] = 0;
  349.     }
  350.     }
  351.  
  352. /* process effects parameters (after delay line initializations, because
  353.     delay line tap ptrs set up here */
  354. reverb->stereoize = stereoize;
  355. reverb->snapToPrime = snapToPrime;
  356. reverb->reverberationTime = reverberationTime;
  357. ProcessParameters(effectType, samplingRate, MAX_CHANNELS, reverb);
  358.  
  359. if (monoProcess == TRUE)
  360.     channelsToProcess = 1;
  361. else
  362.     channelsToProcess = channelCount;
  363.  
  364. /* input audio data, process, output */
  365. while (TRUE)
  366.     {
  367.     ALreadsamps(audioInPort, shortBuffer, BLOCK_SIZE*channelCount);
  368.     if (effectLevel != 0)
  369.     {
  370.     DeinterleaveShortsToNFloats(shortBuffer, systemInputs, BLOCK_SIZE, channelCount); 
  371.     
  372.     if (channelsToProcess == 1)
  373.         {
  374.         SumNFloats(systemInputs, tmpBuffer[0], BLOCK_SIZE, channelCount);
  375.         in = &tmpBuffer[0];
  376.         }
  377.     else
  378.         in = systemInputs;
  379.  
  380.     switch (effectType)
  381.         {
  382.         case SINGLE_ALL_PASS:
  383.         case SINGLE_COMB:
  384.         iiirType = 1;
  385.         if (effectType == SINGLE_ALL_PASS)
  386.             iiirType = 3;
  387.  
  388.         for (i = 0; i < channelsToProcess; i++)
  389.             {
  390.         /* compute comb or all-pass filter */
  391.             IIIRFloats(in[i], reverb->delayOutputs[i][0], BLOCK_SIZE,
  392.                     reverb->delayLines[i][0], reverb->delayLineLengths[i][0], 
  393.                     &reverb->delayLinePtrs[i][0], reverb->delayLoopGains[i][0], iiirType);
  394.             ScaleFloats(reverb->delayOutputs[i][0], reverb->delayOutputs[i][0], BLOCK_SIZE, effectLevel);
  395.             }
  396.  
  397.         /* mix processed signal w/input */
  398.         for (i = 0, delayChannel = 0; i < channelCount; i++)
  399.             {
  400.             ScaleFloats(systemInputs[i], systemInputs[i], BLOCK_SIZE, 1-effectLevel);
  401.             SumFloats(reverb->delayOutputs[delayChannel][0], systemInputs[i], outputs[i], BLOCK_SIZE);
  402.  
  403.             if (channelsToProcess > 1)
  404.             delayChannel++;
  405.             }
  406.         break;
  407.     
  408.         case QUINT_ALL_PASS:
  409.     /* compute 5 series all-pass filters */
  410.         for (i = 0; i < channelsToProcess; i++)
  411.             {
  412.             inPtr = in[i];
  413.             for (j = 0; j < 5; j++)
  414.             {
  415.             IIIRFloats(inPtr, reverb->delayOutputs[i][j], BLOCK_SIZE,
  416.                     reverb->delayLines[i][j], reverb->delayLineLengths[i][j], 
  417.                     &reverb->delayLinePtrs[i][j], reverb->delayLoopGains[i][j], 3);
  418.             inPtr = reverb->delayOutputs[i][j];
  419.             }
  420.  
  421.             ScaleFloats(reverb->delayOutputs[i][j-1], reverb->delayOutputs[i][j-1], BLOCK_SIZE, effectLevel);            
  422.             }
  423.         
  424.         /* mix processed signal w/input */
  425.         for (i = 0, delayChannel = 0; i < channelCount; i++)
  426.             {
  427.             ScaleFloats(systemInputs[i], systemInputs[i], BLOCK_SIZE, 1-effectLevel);
  428.             SumFloats(reverb->delayOutputs[delayChannel][j-1], systemInputs[i], outputs[i], BLOCK_SIZE);
  429.  
  430.             if (channelsToProcess > 1)
  431.             delayChannel++;
  432.             }
  433.         break;
  434.     
  435.         case QUINT_ALL_PASS_STEREO:
  436.     /* compute 5 series all-pass filters */
  437.         for (i = 0; i < channelsToProcess; i++)
  438.             {
  439.             inPtr = in[i];
  440.             for (j = 0; j < 5; j++)
  441.             {
  442.             IIIRFloats(inPtr, reverb->delayOutputs[i][j], BLOCK_SIZE,
  443.                     reverb->delayLines[i][j], reverb->delayLineLengths[i][j], 
  444.                     &reverb->delayLinePtrs[i][j], reverb->delayLoopGains[i][j], 3);
  445.             inPtr = reverb->delayOutputs[i][j];
  446.             }
  447.  
  448.             ScaleFloats(reverb->delayOutputs[i][j-1], reverb->delayOutputs[i][j-1], BLOCK_SIZE, effectLevel);            
  449.             }
  450.         
  451.         /* mix processed signal w/input */
  452.         for (i = 0, delayChannel = 0; i < channelCount; i++)
  453.             {
  454.             ScaleFloats(systemInputs[i], systemInputs[i], BLOCK_SIZE, 1-effectLevel);
  455.             SumFloats(reverb->delayOutputs[delayChannel][j-1], systemInputs[i], outputs[i], BLOCK_SIZE);
  456.  
  457.             if (channelsToProcess > 1)
  458.             delayChannel++;
  459.             }
  460.         break;
  461.  
  462.         case QUAD_COMB_DUAL_ALL_PASS:
  463.         for (i = 0; i < channelsToProcess; i++)
  464.             {
  465.         /* compute 4 parallel comb filters & sum*/
  466.             for (j = 0; j < 4; j++)
  467.             IIIRFloats(in[i], reverb->delayOutputs[i][j], BLOCK_SIZE,
  468.                         reverb->delayLines[i][j], reverb->delayLineLengths[i][j], 
  469.                         &reverb->delayLinePtrs[i][j], reverb->delayLoopGains[i][j], 1);
  470.             SumNFloats(reverb->delayOutputs[i], reverb->delayOutputs[i][j-1], BLOCK_SIZE, j);
  471.  
  472.         /* pass sum thru 2 series all-pass filters */
  473.             for (; j < 6; j++)
  474.             IIIRFloats(reverb->delayOutputs[i][j-1], reverb->delayOutputs[i][j], BLOCK_SIZE,
  475.                     reverb->delayLines[i][j], reverb->delayLineLengths[i][j], 
  476.                     &reverb->delayLinePtrs[i][j], reverb->delayLoopGains[i][j], 3);
  477.  
  478.             ScaleFloats(reverb->delayOutputs[i][j-1], reverb->delayOutputs[i][j-1], BLOCK_SIZE, effectLevel);            
  479.             }
  480.         
  481.         /* mix processed signal w/input */
  482.         for (i = 0, delayChannel = 0; i < channelCount; i++)
  483.             {
  484.             ScaleFloats(systemInputs[i], systemInputs[i], BLOCK_SIZE, 1-effectLevel);
  485.             SumFloats(reverb->delayOutputs[delayChannel][j-1], systemInputs[i], outputs[i], BLOCK_SIZE);
  486.  
  487.             if (channelsToProcess > 1)
  488.             delayChannel++;
  489.             }
  490.         break;
  491.  
  492.         case HEX_COMB_ALL_PASS:
  493.         for (i = 0; i < channelsToProcess; i++)
  494.             {
  495.         /* compute 6 parallel comb filters & sum */
  496.             for (j = 0; j < 6; j++)
  497.             IIIRFloats(in[i], reverb->delayOutputs[i][j], BLOCK_SIZE,
  498.                     reverb->delayLines[i][j], reverb->delayLineLengths[i][j], 
  499.                         &reverb->delayLinePtrs[i][j], reverb->delayLoopGains[i][j], 1);
  500.             SumNFloats(reverb->delayOutputs[i], reverb->delayOutputs[i][j-1], BLOCK_SIZE, j);
  501.  
  502.         /* pass sum thru single all-pass filter */
  503.             IIIRFloats(reverb->delayOutputs[i][j-1], reverb->delayOutputs[i][j], BLOCK_SIZE,
  504.                 reverb->delayLines[i][j], reverb->delayLineLengths[i][j], 
  505.                     &reverb->delayLinePtrs[i][j], reverb->delayLoopGains[i][j], 3);
  506.  
  507.             ScaleFloats(reverb->delayOutputs[i][j], reverb->delayOutputs[i][j], BLOCK_SIZE, effectLevel);            
  508.             }
  509.         
  510.         /* mix processed signal w/input */
  511.         for (i = 0, delayChannel = 0; i < channelCount; i++)
  512.             {
  513.             ScaleFloats(systemInputs[i], systemInputs[i], BLOCK_SIZE, 1-effectLevel);
  514.             SumFloats(reverb->delayOutputs[delayChannel][j], systemInputs[i], outputs[i], BLOCK_SIZE);
  515.  
  516.             if (channelsToProcess > 1)
  517.             delayChannel++;
  518.             }
  519.         break;
  520.     
  521.         case ER7:
  522.         case ER19:
  523.         for (i = 0; i < channelsToProcess; i++)
  524.             {
  525.         /* compute early reflection pattern */
  526.             tapCount = 7;
  527.             if (effectType == ER19)
  528.             tapCount = 19;
  529.             MultiTapDelayFloats(in[i], reverb->delayOutputs[i][0], BLOCK_SIZE, 
  530.                     reverb->delayLines[i][0], reverb->delayLineLengths[i][tapCount-1]+1,
  531.                     reverb->delayLoopGains[i], reverb->delayLinePtrs[i], tapCount);
  532.             }
  533.         
  534.         /* mix processed signal w/input */
  535.         for (i = 0, delayChannel = 0; i < channelCount; i++)
  536.             {
  537.             CopyFloats(reverb->delayOutputs[delayChannel][0], outputs[i], BLOCK_SIZE);
  538.             if (channelsToProcess > 1)
  539.             delayChannel++;
  540.             }
  541.         break;
  542.  
  543.         case HEX_COMB_ALL_PASS_LPF:
  544.         for (i = 0; i < channelsToProcess; i++)
  545.             {
  546.         /* compute 6 parallel comb filters & sum */
  547.             for (j = 0; j < 6; j++)
  548.             IIIRFloatsLPF(in[i], reverb->delayOutputs[i][j], BLOCK_SIZE,
  549.                     reverb->delayLines[i][j], reverb->delayLineLengths[i][j], 
  550.                         &reverb->delayLinePtrs[i][j], reverb->delayLoopGains[i][j], 
  551.                         &reverb->inLineDelays[i][j], reverb->inLineDelayGains[i][j], 1);
  552.             SumNFloats(reverb->delayOutputs[i], reverb->delayOutputs[i][j-1], BLOCK_SIZE, j);
  553.  
  554.         /* pass sum thru single all-pass filter */
  555.             IIIRFloats(reverb->delayOutputs[i][j-1], reverb->delayOutputs[i][j], BLOCK_SIZE,
  556.                 reverb->delayLines[i][j], reverb->delayLineLengths[i][j], 
  557.                     &reverb->delayLinePtrs[i][j], reverb->delayLoopGains[i][j], 3);
  558.  
  559.             ScaleFloats(reverb->delayOutputs[i][j], reverb->delayOutputs[i][j], BLOCK_SIZE, effectLevel);            
  560.             }
  561.         
  562.         /* mix processed signal w/input */
  563.         for (i = 0, delayChannel = 0; i < channelCount; i++)
  564.             {
  565.             ScaleFloats(systemInputs[i], systemInputs[i], BLOCK_SIZE, 1-effectLevel);
  566.             SumFloats(reverb->delayOutputs[delayChannel][j], systemInputs[i], outputs[i], BLOCK_SIZE);
  567.  
  568.             if (channelsToProcess > 1)
  569.             delayChannel++;
  570.             }
  571.         break;
  572.  
  573.         case HEX_COMB_ALL_PASS_LPF_ER7:
  574.         case HEX_COMB_ALL_PASS_LPF_ER19:
  575.         for (i = 0; i < channelsToProcess; i++)
  576.             {
  577.         /* compute early reflection pattern */
  578.             tapCount = 19;
  579.             if (effectType == HEX_COMB_ALL_PASS_LPF_ER7)
  580.             tapCount = 7;
  581.             MultiTapDelayFloats(in[i], reverb->delayOutputs[i][7], BLOCK_SIZE, 
  582.                     reverb->delayLines[i][7], reverb->delayLineLengths[i][7+tapCount-1]+1,
  583.                     &reverb->delayLoopGains[i][7], &reverb->delayLinePtrs[i][7], tapCount);
  584.  
  585.         /* compute 6 parallel comb filters & sum */
  586.             for (j = 0; j < 6; j++)
  587.             IIIRFloatsLPF(reverb->delayOutputs[i][7], reverb->delayOutputs[i][j], BLOCK_SIZE,
  588.                     reverb->delayLines[i][j], reverb->delayLineLengths[i][j], 
  589.                         &reverb->delayLinePtrs[i][j], reverb->delayLoopGains[i][j], 
  590.                         &reverb->inLineDelays[i][j], reverb->inLineDelayGains[i][j], 1);
  591.             SumNFloats(reverb->delayOutputs[i], reverb->delayOutputs[i][5], BLOCK_SIZE, 6);
  592.  
  593.         /* pass sum thru single all-pass filter */
  594.             IIIRFloats(reverb->delayOutputs[i][5], reverb->delayOutputs[i][6], BLOCK_SIZE,
  595.                 reverb->delayLines[i][6], reverb->delayLineLengths[i][6], 
  596.                     &reverb->delayLinePtrs[i][6], reverb->delayLoopGains[i][6], 3);
  597.  
  598. #define LATE_DELAY  7+tapCount
  599.         /* delay late reflections */
  600.             DelayFloats(reverb->delayOutputs[i][6], reverb->delayOutputs[i][LATE_DELAY], BLOCK_SIZE, 
  601.                     reverb->delayLines[i][LATE_DELAY], reverb->delayLineLengths[i][LATE_DELAY], 
  602.                     &reverb->delayLinePtrs[i][LATE_DELAY], FALSE);
  603.  
  604.         /* add delayed late reflections to early reflections */
  605.             SumFloats(reverb->delayOutputs[i][LATE_DELAY], reverb->delayOutputs[i][7], 
  606.                 reverb->delayOutputs[i][LATE_DELAY], BLOCK_SIZE);
  607.             ScaleFloats(reverb->delayOutputs[i][LATE_DELAY], reverb->delayOutputs[i][LATE_DELAY], BLOCK_SIZE, effectLevel);            
  608.             }
  609.         
  610.         /* mix processed signal w/input */
  611.         for (i = 0, delayChannel = 0; i < channelCount; i++)
  612.             {
  613.             ScaleFloats(systemInputs[i], systemInputs[i], BLOCK_SIZE, 1-effectLevel);
  614.             SumFloats(reverb->delayOutputs[delayChannel][LATE_DELAY], systemInputs[i], outputs[i], BLOCK_SIZE);
  615.             if (channelsToProcess > 1)
  616.             delayChannel++;
  617.             }
  618.  
  619.         break;
  620.  
  621.         default:
  622.         break;
  623.         }
  624.  
  625.     InterleaveNFloatsToShorts(outputs, shortBuffer, BLOCK_SIZE, channelCount, TRUE);
  626.     }
  627.  
  628.     ALwritesamps(audioOutPort, shortBuffer, BLOCK_SIZE*channelCount);
  629.     }
  630.  
  631. /* release audio ports */
  632. ALcloseport(audioInPort);
  633. ALcloseport(audioOutPort);
  634. }    /* ------------------ end main() --------------- */
  635.  
  636. /*
  637.  * GetHardwareInputRate:    acquire audio hardware input sampling rate
  638.  * -------------------------------------------------------------------- */
  639.     long 
  640. GetHardwareInputRate()
  641. {
  642. long buffer[6];
  643.  
  644. /* acquire state variables of audio hardware */
  645. buffer[0] = AL_INPUT_RATE;
  646. buffer[2] = AL_INPUT_SOURCE;
  647. buffer[4] = AL_DIGITAL_INPUT_RATE;
  648. ALgetparams(AL_DEFAULT_DEVICE, buffer, 6);
  649.  
  650. /* for input sources microphone or line and input rate not AES word clock */
  651. if    ((buffer[3] != AL_INPUT_DIGITAL)&&(buffer[1] > 0))
  652.     return (buffer[1]);
  653.  
  654. /* for input rate AES word clock and machine has ability to read digital input sampling rate, return 
  655.     AL_DIGITAL_INPUT_RATE */
  656. else if (ALgetdefault(AL_DEFAULT_DEVICE, AL_DIGITAL_INPUT_RATE) >= 0) 
  657.     return (buffer[5]);
  658.  
  659. /* could not determine sampling rate */
  660. else
  661.     return (AL_RATE_UNDEFINED);
  662. }    /* ------------------------- end GetHardwareInputRate() -------------------------- */
  663.  
  664. /* **********************************************************************
  665.  * ParseCommandLine:    parse input command line for user options  
  666.  * ********************************************************************** */
  667.     static void
  668. ParseCommandLine(int argc, char **argv, char fileName[][200])
  669. {
  670. int        i;
  671.  
  672. for (i = 1; i < argc; i++)
  673.     {
  674. /* check for -h or -help */
  675.     if ((!strcmp(argv[i], "-h"))||(!strcmp(argv[i], "-help")))
  676.     {
  677.     fprintf(stderr, "Usage: %s\n", applicationUsage);
  678.     exit(1);
  679.     }
  680.  
  681. /* check for -verbose */
  682.     else if ((!strcmp(argv[i], "-v"))||(!strcmp(argv[i], "-verbose")))
  683.     {
  684.     verboseOperation = TRUE;
  685.     }
  686.  
  687. /* check for -noMonoProcess */
  688.     else if (!strcmp(argv[i], "-noMonoProcess"))
  689.     {
  690.     monoProcess = FALSE;
  691.     }
  692.  
  693. /* check for -noPrime */
  694.     else if (!strcmp(argv[i], "-noPrime"))
  695.     {
  696.     snapToPrime = FALSE;
  697.     }
  698.  
  699. /* check for reverbTime */
  700.     else if (!strcmp(argv[i], "-reverbTime"))
  701.     {
  702.     i++;
  703.     if (i < argc)
  704.         {
  705.         reverberationTime = atof(argv[i]);
  706.         if (reverberationTime < 0)
  707.         {
  708.         fprintf(stderr, "Bogus %s: %s\n", argv[i-1], argv[i]);
  709.         exit(1);
  710.         }
  711.         }
  712.     else
  713.         {
  714.         fprintf(stderr, "Usage: %s\n", applicationUsage);
  715.         exit(1);
  716.         }
  717.     }
  718.  
  719. /* check for effectLevel */
  720.     else if (!strcmp(argv[i], "-level"))
  721.     {
  722.     i++;
  723.     if (i < argc)
  724.         {
  725.         effectLevel = atof(argv[i]);
  726.         if (effectLevel < 0)
  727.         {
  728.         fprintf(stderr, "Bogus %s: %s\n", argv[i-1], argv[i]);
  729.         exit(1);
  730.         }
  731.         }
  732.     else
  733.         {
  734.         fprintf(stderr, "Usage: %s\n", applicationUsage);
  735.         exit(1);
  736.         }
  737.     }
  738.  
  739. /* check for effect type */
  740.     else if (!strcmp(argv[i], "-type"))
  741.     {
  742.     i++;
  743.     if (i < argc)
  744.         {
  745.         if        (!strcmp(argv[i], "schroeder1"))
  746.         effectType = QUINT_ALL_PASS;
  747.         else if (!strcmp(argv[i], "schroeder2"))
  748.         effectType = QUAD_COMB_DUAL_ALL_PASS;
  749.  
  750.         else if (!strcmp(argv[i], "chamberlin"))
  751.         effectType = QUINT_ALL_PASS_STEREO;
  752.  
  753.         else if (!strcmp(argv[i], "moorer"))
  754.         effectType = HEX_COMB_ALL_PASS;
  755.         else if (!strcmp(argv[i], "moorer2"))
  756.         effectType = HEX_COMB_ALL_PASS_LPF;
  757.  
  758.         else if (!strcmp(argv[i], "er7"))
  759.         effectType = ER7;
  760.         else if (!strcmp(argv[i], "er19"))
  761.         effectType = ER19;
  762.         else if (!strcmp(argv[i], "moorer2er7"))
  763.         effectType = HEX_COMB_ALL_PASS_LPF_ER7;
  764.         else if (!strcmp(argv[i], "moorer2er19"))
  765.         effectType = HEX_COMB_ALL_PASS_LPF_ER19;
  766.  
  767.         else if (!strcmp(argv[i], "comb"))
  768.         effectType = SINGLE_COMB;
  769.         else if (!strcmp(argv[i], "allpass"))
  770.         effectType = SINGLE_ALL_PASS;
  771.  
  772.         else
  773.         {
  774.         fprintf(stderr, "Bogus %s: %s\n", argv[i-1], argv[i]);
  775.         exit(1);
  776.         }
  777.         }
  778.     else
  779.         {
  780.         fprintf(stderr, "Usage: %s\n", applicationUsage);
  781.         exit(1);
  782.         }
  783.     }
  784.  
  785.     /* oh, error */
  786.     else
  787.     {
  788.     fprintf(stderr, "Usage: %s\n", applicationUsage);
  789.     exit(1);
  790.     }
  791.     }
  792. }    /* ----------------------- end ParseCommandLine() --------------------- */
  793.  
  794. /* **********************************************************************
  795.  * ProcessParameters:    parse input command line for user options  
  796.  * ********************************************************************** */
  797.     static void
  798. ProcessParameters(char type, double samplingRate, int channels, Reverb *reverb)
  799. /*
  800.     type            algorithm type
  801.     samplingRate        sampling rate
  802.     channels            # processing channels
  803.     reverb            ptr to parameter structure
  804. */
  805. {
  806. int        i, j, k;
  807. float        hiCoeff, loCoeff;
  808. double        *times;
  809. float        *gains, *inLineGains;
  810. int        *lengths, *delayLinePtrs;
  811. int        channel;
  812. int        value;
  813. double        reverberationTime;
  814. int        tapCount;
  815. int        maxLength;
  816. double        times_er[19], gains_er[19];
  817. int        earlyToLateReflectionAlignmentLength;
  818.  
  819. reverberationTime = reverb->reverberationTime;
  820.  
  821. /* initialize early reflections parameters */
  822. if    ((type == ER7)||(type == HEX_COMB_ALL_PASS_LPF_ER7))
  823.     {
  824.     times_er[0] = 0;        gains_er[0] = 1;
  825.     times_er[1] = 0.0199;   gains_er[1] = 1.02;
  826.     times_er[2] = 0.0354;   gains_er[2] = 0.818;
  827.     times_er[3] = 0.0389;   gains_er[3] = 0.635;
  828.     times_er[4] = 0.0414;   gains_er[4] = 0.719;
  829.     times_er[5] = 0.0699;   gains_er[5] = 0.267;
  830.     times_er[6] = 0.0796;   gains_er[6] = 0.242;
  831.     }
  832. else if ((type == ER19)||(type == HEX_COMB_ALL_PASS_LPF_ER19))
  833.     {
  834.     times_er[0]  = 0;         gains_er[0]  = 1;
  835.     times_er[1]  = 0.0043;   gains_er[1]  = 0.841;
  836.     times_er[2]  = 0.0215;   gains_er[2]  = 0.504;
  837.     times_er[3]  = 0.0225;   gains_er[3]  = 0.491;
  838.     times_er[4]  = 0.0268;   gains_er[4]  = 0.379;
  839.     times_er[5]  = 0.0270;   gains_er[5]  = 0.380;
  840.     times_er[6]  = 0.0298;   gains_er[6]  = 0.346;
  841.     times_er[7]  = 0.0458;   gains_er[7]  = 0.289;
  842.     times_er[8]  = 0.0485;   gains_er[8]  = 0.272;
  843.     times_er[9]  = 0.0572;   gains_er[9]  = 0.192;
  844.     times_er[10] = 0.0587;   gains_er[10] = 0.193;
  845.     times_er[11] = 0.0595;   gains_er[11] = 0.217;
  846.     times_er[12] = 0.0612;   gains_er[12] = 0.181;
  847.     times_er[13] = 0.0707;   gains_er[13] = 0.180;
  848.     times_er[14] = 0.0708;   gains_er[14] = 0.181;
  849.     times_er[15] = 0.0726;   gains_er[15] = 0.176;
  850.     times_er[16] = 0.0741;   gains_er[16] = 0.142;
  851.     times_er[17] = 0.0753;   gains_er[17] = 0.167;
  852.     times_er[18] = 0.0797;   gains_er[18] = 0.134;
  853.     }
  854.  
  855. /* 
  856.  * initialize reverberation parameters according to algorithm
  857.  */
  858. switch (type)
  859.     {
  860.     case SINGLE_COMB:
  861.     case SINGLE_ALL_PASS:
  862.     for (channel = 0; channel < channels; channel++)  
  863.         {      
  864.         times = reverb->delayLineTimes[channel];
  865.         gains = reverb->delayLoopGains[channel];
  866.         lengths = reverb->delayLineLengths[channel];
  867.  
  868.         times[0] = 0.100;
  869.         gains[0] = pow(10, -3*times[0]/reverberationTime);
  870.  
  871.     /* place opposite gain polarities in adjacent channels for stereo effect */
  872.         if (((i&1) == 0)&&(reverb->stereoize == TRUE))
  873.         gains[0] = -gains[0];
  874.  
  875.         ConvertDelayLengths(lengths, times, samplingRate, 1, reverb->snapToPrime);
  876.         }
  877.     break;
  878.  
  879.     case QUINT_ALL_PASS:
  880.     for (channel = 0; channel < channels; channel++)  
  881.         {      
  882.         times = reverb->delayLineTimes[channel];
  883.         gains = reverb->delayLoopGains[channel];
  884.         lengths = reverb->delayLineLengths[channel];
  885.  
  886.         times[0] = 0.100;
  887.         times[1] = 0.068;
  888.         times[2] = 0.060;
  889.         times[3] = 0.0197;
  890.         times[4] = 0.00585;
  891.         
  892.         gains[0] = 0.7;
  893.         gains[1] = -0.7;
  894.         gains[2] = 0.7;
  895.         gains[3] = 0.7;
  896.         gains[4] = 0.7;
  897.     
  898.         ConvertDelayLengths(lengths, times, samplingRate, 5, reverb->snapToPrime);
  899.         }
  900.     break;
  901.  
  902.     case QUINT_ALL_PASS_STEREO:
  903.     for (channel = 0; channel < channels; channel++)  
  904.         {      
  905.         times = reverb->delayLineTimes[channel];
  906.         gains = reverb->delayLoopGains[channel];
  907.         lengths = reverb->delayLineLengths[channel];
  908.  
  909.         times[0] = 0.0496;
  910.         times[1] = 0.03465;
  911.         times[2] = 0.02418;
  912.  
  913.         gains[0] = 0.75;
  914.         gains[1] = 0.72;
  915.         gains[2] = 0.691;
  916.  
  917.     /* use Lauridsen effect for stereoization */
  918.         if ((i&1) == 0)
  919.         {
  920.         times[3] = 0.01785;
  921.         times[4] = 0.01098;
  922.         gains[3] = 0.649;
  923.         gains[4] = 0.662;
  924.         }
  925.         else
  926.         {
  927.         times[3] = 0.01801;
  928.         times[4] = 0.01082;
  929.         gains[3] = 0.646;
  930.         gains[4] = 0.666;
  931.         }
  932.         
  933.         ConvertDelayLengths(lengths, times, samplingRate, 5, reverb->snapToPrime);
  934.         }
  935.  
  936.     break;
  937.  
  938.     case QUAD_COMB_DUAL_ALL_PASS:
  939.     for (channel = 0; channel < channels; channel++)  
  940.         {      
  941.         times = reverb->delayLineTimes[channel];
  942.         gains = reverb->delayLoopGains[channel];
  943.         lengths = reverb->delayLineLengths[channel];
  944.  
  945.     /* comb filters (shortest = predelay: time between direct sound and reverb) */
  946.         times[0] = 0.045;    
  947.         times[1] = 0.040;
  948.         times[2] = 0.035;
  949.         times[3] = 0.030;
  950.     
  951.     /* comb filters: gain = 10^(-3*delayTime/reverberationTime) */
  952.         for (j = 0; j < 4; j++)
  953.             gains[j] = pow(10, -3*times[j]/reverberationTime);
  954.     
  955.     /* all-pass filters */
  956.         times[4] = 0.005;
  957.         gains[4] = 0.7;
  958.         times[5] = 0.0017;
  959.         gains[5] = 0.7;
  960.         
  961.         ConvertDelayLengths(lengths, times, samplingRate, 6, reverb->snapToPrime);
  962.         }
  963.     break;
  964.  
  965.     case HEX_COMB_ALL_PASS:
  966.     case HEX_COMB_ALL_PASS_LPF:
  967.     for (channel = 0; channel < channels; channel++)  
  968.         {      
  969.         times = reverb->delayLineTimes[channel];
  970.         gains = reverb->delayLoopGains[channel];
  971.         lengths = reverb->delayLineLengths[channel];
  972.         inLineGains = reverb->inLineDelayGains[channel];
  973.  
  974.     /* comb filters (shortest = predelay: time between direct sound and reverb) */
  975.         times[0] = 0.050;    
  976.         times[1] = 0.056;
  977.         times[2] = 0.061;
  978.         times[3] = 0.068;
  979.         times[4] = 0.072;
  980.         times[5] = 0.078;
  981.         
  982.     /* comb filters: gain = 10^(-3*delayTime/reverberationTime) */
  983.     /* gain = 10^(-3*delayTime/reverberationTime) but moorer sets all to same gain 
  984.         with time 0.050 seconds */
  985.         for (j = 0; j < 6; j++)
  986.         gains[j] = pow(10, -3*0.050/reverberationTime);
  987.     
  988.     /* all-pass filters */
  989.         times[6] = 0.006;
  990.         gains[6] = 0.7;
  991.         
  992.         if (type == HEX_COMB_ALL_PASS_LPF)
  993.         {
  994.         /* in-loop low pass filter coefficients */
  995.         /* linearly interpolate new point from 2 tables at 25000 & 50000 Hz 
  996.             y = y0 + (x-x0)*(y1-y0)/(x1-x0) */
  997.         for (i = 0; i < 6; i++)
  998.             {
  999.             hiCoeff = highFrequencyDampingCoeffs_50000Hz[i];
  1000.             loCoeff = highFrequencyDampingCoeffs_25000Hz[i];
  1001.             inLineGains[i] = loCoeff +  
  1002.             (samplingRate - 25000)*(hiCoeff - loCoeff)/(50000 - 25000);
  1003.  
  1004.         /* force to stability condition:  g2/(1-g1) < 1 --> g2 = g(1 - g1) */
  1005.         /* gains[j] = (1 - inLineGains[j])*pow(10, -3*times[j]/reverberationTime);*/
  1006.             gains[i] *= (1 - inLineGains[i]);
  1007.             }    
  1008.         }        
  1009.      
  1010.         ConvertDelayLengths(lengths, times, samplingRate, 7, reverb->snapToPrime);
  1011.         }
  1012.     break;
  1013.    
  1014.     case ER7:
  1015.     case ER19:
  1016.     for (channel = 0; channel < channels; channel++)  
  1017.         {      
  1018.         times = reverb->delayLineTimes[channel];
  1019.         gains = reverb->delayLoopGains[channel];
  1020.         lengths = reverb->delayLineLengths[channel];
  1021.         delayLinePtrs = reverb->delayLinePtrs[channel];
  1022.  
  1023.         if (type == ER7)
  1024.         tapCount = 7;
  1025.         else
  1026.         tapCount = 19;
  1027.         for (i = 0; i < tapCount; i++)
  1028.         {
  1029.         times[i] = times_er[i];
  1030.         gains[i] = gains_er[i];
  1031.         }        
  1032.         ConvertDelayLengths(lengths, times, samplingRate, tapCount, FALSE);
  1033.  
  1034.     /* initialize tap ptrs.  Longest tap assumed to be last tap */      
  1035.         delayLinePtrs[i] = 0;  
  1036.         maxLength = lengths[tapCount-1];
  1037.         for (i = 1; i < tapCount; i++)
  1038.         delayLinePtrs[i] = 1 + maxLength - lengths[i];
  1039.         }
  1040.     break;
  1041.  
  1042.     case HEX_COMB_ALL_PASS_LPF_ER7:
  1043.     case HEX_COMB_ALL_PASS_LPF_ER19:
  1044.     for (channel = 0; channel < channels; channel++)  
  1045.         {      
  1046.         times = reverb->delayLineTimes[channel];
  1047.         gains = reverb->delayLoopGains[channel];
  1048.         lengths = reverb->delayLineLengths[channel];
  1049.         inLineGains = reverb->inLineDelayGains[channel];
  1050.         delayLinePtrs = reverb->delayLinePtrs[channel];
  1051.  
  1052.     /* 
  1053.      * LATE REFLECTION PATTERN 
  1054.      */
  1055.     /* comb filters (shortest = predelay: time between direct sound and reverb) */
  1056.         times[0] = 0.050;    
  1057.         times[1] = 0.056;
  1058.         times[2] = 0.061;
  1059.         times[3] = 0.068;
  1060.         times[4] = 0.072;
  1061.         times[5] = 0.078;
  1062.         
  1063.     /* all-pass filter */
  1064.         times[6] = 0.006;
  1065.         gains[6] = 0.7;
  1066.         
  1067.         for (i = 0; i < 6; i++)
  1068.         {
  1069.         /* comb filters: gain = 10^(-3*delayTime/reverberationTime) */
  1070.         /* gain = 10^(-3*delayTime/reverberationTime) (Moorer sets all to same gain) */
  1071.         /* gain = (1 - inLineGain)*pow(10, -3*time/reverberationTime);*/
  1072.         gains[i] = pow(10, -3*0.05/reverberationTime);
  1073.  
  1074.         /* in-loop low pass filter coefficients */
  1075.         /* linearly interpolate new point from 25000 & 50000 Hz tables */
  1076.         /*y = y0 + (x-x0)*(y1-y0)/(x1-x0) */
  1077.         hiCoeff = highFrequencyDampingCoeffs_50000Hz[i];
  1078.         loCoeff = highFrequencyDampingCoeffs_25000Hz[i];
  1079.         inLineGains[i] = loCoeff +  
  1080.             (samplingRate - 25000)*(hiCoeff - loCoeff)/(50000 - 25000);
  1081.  
  1082.         /* force to stability condition:  g2/(1-g1) < 1 --> g2 = g(1 - g1) */
  1083.         gains[i] *= (1 - inLineGains[i]);
  1084.         }
  1085.      
  1086.         ConvertDelayLengths(lengths, times, samplingRate, 7, reverb->snapToPrime);
  1087.  
  1088.     /* 
  1089.      * EARLY REFLECTION PATTERN
  1090.      */
  1091.         tapCount = 7;
  1092.         if (type == HEX_COMB_ALL_PASS_LPF_ER19)
  1093.         tapCount = 19;
  1094.         times += 7;
  1095.         gains += 7;
  1096.         lengths += 7;
  1097.         delayLinePtrs += 7;
  1098.  
  1099.         for (i = 0; i < tapCount; i++)
  1100.         {
  1101.         times[i] = times_er[i];
  1102.         gains[i] = gains_er[i];
  1103.         }      
  1104.  
  1105.         ConvertDelayLengths(lengths, times, samplingRate, tapCount, FALSE);
  1106.  
  1107.     /*  align last tap of ER generator lands close to first echo from 
  1108.         late reflection by adding delay to late reflections
  1109.         alignment delay = ER longest delay - LR shortest comb delay 
  1110.  
  1111.         assume last ER delay is longer than shortest LR comb delay, as is
  1112.         case for 7 and 19 tap ER generators       
  1113.         */
  1114.         lengths[tapCount] = lengths[tapCount-1] - lengths[-7];
  1115.  
  1116.     /* initialize tap ptrs.  Longest tap assumed to be last tap */      
  1117.         delayLinePtrs[i] = 0;  
  1118.         maxLength = lengths[tapCount-1];
  1119.         for (i = 0; i < tapCount; i++)
  1120.         delayLinePtrs[i] = 1 + maxLength - lengths[i];
  1121.         }
  1122.     break;
  1123.  
  1124.     default:
  1125.     printf("ProcessParameters(): bogus type: %d\n", type);
  1126.     break;
  1127.     }
  1128.  
  1129.  
  1130.  
  1131. times = reverb->delayLineTimes[0];
  1132. gains = reverb->delayLoopGains[0];
  1133. lengths = reverb->delayLineLengths[0];
  1134. inLineGains = reverb->inLineDelayGains[0];
  1135. delayLinePtrs = reverb->delayLinePtrs[0];
  1136.  
  1137. #ifdef SAFE
  1138. for (i = 0; i < 30; i++)
  1139.     printf("%d: t=%f, g=%f, ig=%f, l=%d\n", i, times[i], gains[i], inLineGains[i], lengths[i]);
  1140. #endif
  1141. } /* ----------------------- end ProcessParameters() --------------------- */
  1142.  
  1143. /* **********************************************************************
  1144.  * IIIRFloatsLPF:   interpolated IIR (comb/all-pass) w/low pass filter
  1145.  * ********************************************************************** */
  1146.     void
  1147. IIIRFloatsLPF(float *input, float *output, int length, 
  1148.             float *delayLine, int delayLength, int *delayLineIndex, float loopGain, 
  1149.             float *inLoopDelay, float inLoopGain, 
  1150.                 int topology)
  1151. /*  input        ptr to input buffer
  1152.     output        ptr to output buffer
  1153.     length        length of input buffer 
  1154.     delayLine        ptr to delay line
  1155.     delayLength        length of delay buffer
  1156.     delayLineIndex    ptr to index of next read/write position in delay buffer.
  1157.     loopGain        feedback gain
  1158.     inLoopDelay        in-loop IIR filter delay element
  1159.     inLoopGain        in-loop IIR filter coefficient
  1160.     topology        0:  comb filter,    input not delayed
  1161.             1:  comb filter,    input delayed
  1162.             2:  comb filter,    input delayed, scaled
  1163.             3:  all-pass filter, 
  1164.             4:  all-pass filter, one multiply 
  1165. */
  1166. {
  1167. int    i;
  1168. int    delayIndex;
  1169. float    y, sum;
  1170. float    z;
  1171.  
  1172. delayIndex = *delayLineIndex;
  1173. z = *inLoopDelay;
  1174.  
  1175. switch (topology)
  1176.     {
  1177. /* comb filter: input delayed, IIR low-pass filter in loop */
  1178. /*         -m            -m
  1179.             z           z
  1180.     H(z) = --------------- = ------------
  1181.                 -m        -m
  1182.         1 - g*T(z)*z     1 - g*z    
  1183.                 ---------
  1184.                     -1
  1185.                 1 - g1*z
  1186.  
  1187.     where T(z) = 1/(1 - g1*(z^-1))
  1188. */
  1189.     case 1:
  1190.     for (i = 0; i < length; i++)
  1191.         {
  1192.         y = delayLine[delayIndex];
  1193.         output[i] = y;
  1194.  
  1195.     /* update delay elements */
  1196.         y += z*inLoopGain;
  1197.         delayLine[delayIndex++] = input[i] + y*loopGain;
  1198.         if (delayIndex >= delayLength)
  1199.         delayIndex = 0;
  1200.         z = y;
  1201.         }
  1202.     break;
  1203.  
  1204.     default:
  1205.     fprintf(stderr, "IIIRFloatsLPF(): bogus topology: %d\n", topology);
  1206.     break;
  1207.     }
  1208.  
  1209. *inLoopDelay = z;
  1210. *delayLineIndex = delayIndex;
  1211. }    /* ----------------------- end IIIRFloatsLPF() --------------------- */
  1212.  
  1213. /* **********************************************************************
  1214.  * MultiTapDelayFloats:   interpolated IIR (comb/all-pass) w/low pass filter
  1215.  * ********************************************************************** */
  1216.     static void
  1217. MultiTapDelayFloats(float *input, float *output, int length, 
  1218.             float *delayLine, int delayLength,
  1219.                 float *tapGains, int *tapIndices, int tapCount)
  1220. /*  input        ptr to input buffer
  1221.     output        ptr to output buffer
  1222.     length        length of input buffer 
  1223.     delayLine        ptr to delay line
  1224.     delayLength        length of longest tap
  1225.  
  1226.     tapGains        ptr to delay line tap gains
  1227.     tapIndices        ptr to delay line tap indices
  1228.     tapCount        # taps
  1229. */
  1230. {
  1231. int        i, tap, tapIndex;
  1232. float        sum;
  1233.  
  1234. /* sum taps, advance and bind tap ptrs.  In anticipation of conditional move 
  1235. instruction, change code to decrement counter and wrap at 0 rather than 
  1236. increment and wrap at 'delayLength' */
  1237. for (i = 0; i < length; i++)
  1238.     {
  1239.     delayLine[tapIndices[0]] = input[i];
  1240.     for (tap = 0, sum = 0; tap < tapCount; tap++)
  1241.     {
  1242.     sum += tapGains[tap]*delayLine[tapIndices[tap]++];
  1243.     if (tapIndices[tap] >= delayLength)
  1244.         tapIndices[tap] = 0;
  1245.     }
  1246.     output[i] = sum;
  1247.     }
  1248. }    /* ----------------------- end MultiTapDelayFloats() --------------------- */
  1249.  
  1250. /* **********************************************************************
  1251.  * ConvertDelayLengths:        convert delay times to delay lengths
  1252.  * ********************************************************************** */
  1253.     static void
  1254. ConvertDelayLengths(int *lengths, double *times, double samplingRate,
  1255.             int count, char snapToPrime)
  1256. /*
  1257. */
  1258. {
  1259. int    i;
  1260.  
  1261. for (i = 0; i < count; i++)
  1262.     {
  1263.     lengths[i] = SecondsToSamples(times[i], samplingRate);
  1264.     if (snapToPrime == TRUE)
  1265.     {
  1266.     lengths[i] = NearestPrime(lengths[i], 2, 10000);
  1267.     }
  1268.     }
  1269. }    /* ----------------------- end ConvertDelayLengths() --------------------- */
  1270.  
  1271. /**************************************************************** 
  1272.  * SetFloats   fill buffer w/'value' FLOATs
  1273.  *****************************************************************/
  1274.     void 
  1275. SetFloats(float *in, int length, float value)
  1276. /* in        ptr to input data
  1277.    length   length of data 
  1278.    value    value to fill buffer
  1279. */
  1280. {
  1281. register int     i;            
  1282. for (i = 0; i < length; i++) 
  1283.     in[i] = value;
  1284. }    /* ------------------ end SetFloats() --------------- */
  1285.  
  1286. /* **************************************************************
  1287.  * DeinterleaveShortsToNFloats    deinterleave to N buffers
  1288.  *                
  1289.  *            The first word of input buffer at 'input' 
  1290.  *            will be first word in output buffer 'outputLeft.'
  1291.  *****************************************************************/
  1292.     void 
  1293. DeinterleaveShortsToNFloats(short *input, float *outputs[], int outputLength,
  1294.                 int interleaveFactor)
  1295. /* input        ptr to input buffer
  1296.    outputs        ptr to array of output buffers
  1297.    outputLength        length of output buffer 
  1298.    interleaveFactor    # output buffers (interleave factor)
  1299. */
  1300. {
  1301. register int    i, j;
  1302. register short    *in;
  1303.  
  1304. /* deinterleave of N-sample frames */ 
  1305. in = input;
  1306. for (i = 0; i < outputLength; i++) 
  1307.     {
  1308.     for (j = 0; j < interleaveFactor; j++) 
  1309.     outputs[j][i] = (float) in[j];
  1310.     in += interleaveFactor;
  1311.     }
  1312. }    /* ------------------ end DeinterleaveShortsToNFloats() --------------- */
  1313.  
  1314. /****************************************************************
  1315.  * SumNFloats        sum N buffers 
  1316.  *****************************************************************/
  1317.     void 
  1318. SumNFloats(float *inputs[], float *output, int length, int inputChannels)
  1319. /*  inputs        ptr to array of input buffers
  1320.     output        ptr to output buffer
  1321.     length        length of buffers
  1322.     inputChannels    #channels to be summed 
  1323. */
  1324. {
  1325. register int        i, j;
  1326. register float        sum;
  1327.  
  1328. for (i = 0; i < length; i++) 
  1329.     {
  1330.     sum = 0;
  1331.     for (j = 0; j < inputChannels; j++) 
  1332.     sum += inputs[j][i];
  1333.     output[i] = sum;
  1334.     }
  1335. }    /* ------------------ end SumNFloats() --------------- */
  1336.  
  1337. /* **********************************************************************
  1338.  * IIIRFloats:   interpolated IIR (comb/all-pass)
  1339.  * ********************************************************************** */
  1340.     void
  1341. IIIRFloats(float *in, float *out, int length, 
  1342.             float *delayLine, int delayLength, int *delayLineIndex, 
  1343.             float loopGain, int topology)
  1344. /*  in            ptr to input buffer
  1345.     out            ptr to output buffer
  1346.     length        length of input buffer 
  1347.     delayLine        ptr to delay line
  1348.     delayLength        length of delay buffer
  1349.     delayLineIndex    ptr to index of next read/write position in delay buffer.
  1350.     loopGain        feedback gain
  1351.     topology        0:  comb filter,    input not delayed
  1352.             1:  comb filter,    input delayed
  1353.             2:  comb filter,    input delayed, scaled
  1354.             3:  all-pass filter, 
  1355.             4:  all-pass filter, one multiply 
  1356. */
  1357. {
  1358. int    i;
  1359. int    delayIndex;
  1360. float    y, sum;
  1361.  
  1362. delayIndex = *delayLineIndex;
  1363.  
  1364. /*
  1365. comb filter magnitude response:
  1366.  
  1367.     gain at maxima = 1/(1-loopGain)
  1368.     gain at minima = 1/(1+loopGain)
  1369.  
  1370.     comb period = 1/delayLength
  1371. */
  1372. switch (topology)
  1373.     {
  1374. /* comb filter: input not delayed */
  1375. /*         
  1376.             1
  1377.     H(z) = ----------
  1378.            -m
  1379.         1 - g*z       
  1380. */
  1381.     case 0:
  1382.     for (i = 0; i < length; i++)
  1383.         {
  1384.         y = in[i] + delayLine[delayIndex]*loopGain;
  1385.         out[i] = y;
  1386.     
  1387.         delayLine[delayIndex++] = y;
  1388.         if (delayIndex >= delayLength)
  1389.         delayIndex = 0;
  1390.         }
  1391.     break;
  1392.  
  1393. /* comb filter: input delayed */
  1394. /*         -m
  1395.             z
  1396.     H(z) = ----------
  1397.            -m
  1398.         1 - g*z       
  1399. */
  1400.     case 1:
  1401.     for (i = 0; i < length; i++)
  1402.         {
  1403.         y = delayLine[delayIndex];
  1404.         out[i] = y;
  1405.         delayLine[delayIndex++] = in[i] + y*loopGain;
  1406.         if (delayIndex >= delayLength)
  1407.         delayIndex = 0;
  1408.         }
  1409.     break;
  1410.  
  1411. /* comb filter: input delayed, scaled */
  1412.     case 2:
  1413.     for (i = 0; i < length; i++)
  1414.         {
  1415.         y = delayLine[delayIndex]*loopGain;
  1416.         out[i] = y;
  1417.         delayLine[delayIndex++] = in[i] + y;
  1418.         if (delayIndex >= delayLength)
  1419.         delayIndex = 0;
  1420.         }
  1421.     break;
  1422.  
  1423. /* all-pass filter */
  1424. /*         -m
  1425.         g +    z
  1426.     H(z) = ----------
  1427.            -m
  1428.         1 + g*z       
  1429. */
  1430.     case 3:
  1431.     for (i = 0; i < length; i++)
  1432.         {
  1433.         y = delayLine[delayIndex];
  1434.         out[i] = y - in[i]*loopGain;
  1435.         delayLine[delayIndex++] = y*loopGain + in[i];
  1436.         if (delayIndex >= delayLength)
  1437.         delayIndex = 0;
  1438.         }
  1439.     break;
  1440.  
  1441. /* one multiply all-pass filter */
  1442. /*         -m
  1443.         g +    z
  1444.     H(z) = ----------
  1445.            -m
  1446.         1 + g*z       
  1447. */
  1448.     case 4:
  1449.     for (i = 0; i < length; i++)
  1450.         {
  1451.         y = delayLine[delayIndex];
  1452.         sum = (in[i] - y)*loopGain;
  1453.  
  1454.         out[i] = sum + y;
  1455.         delayLine[delayIndex++] = sum + in[i];
  1456.         if (delayIndex >= delayLength)
  1457.         delayIndex = 0;
  1458.         }
  1459.     break;
  1460.     default:
  1461.     fprintf(stderr, "IIIRFloats(): bogus topology: %d\n", topology);
  1462.     break;
  1463.     }
  1464. *delayLineIndex = delayIndex;
  1465. }    /* ----------------------- end IIIRFloats() --------------------- */
  1466.  
  1467. /**************************************************************** 
  1468.  * ScaleFloats    scale FLOATs from in buffer to out buffer
  1469.  *                (ok for "in place" scale )
  1470.  *****************************************************************/
  1471.     void 
  1472. ScaleFloats(float *in, float *out, int length, float scaleFactor)
  1473. /* in        ptr to input
  1474.    out        ptr to output
  1475.    length    length of data 
  1476.    scaleFactor    input file to be scaled w/this factor */
  1477. {
  1478. int     i;            
  1479.  
  1480. for (i = 0; i < length; i++)
  1481.     out[i] = in[i]*scaleFactor;
  1482. }    /* ------------------ end ScaleFloats() --------------- */
  1483.  
  1484. /****************************************************************
  1485.  * SumFloats        sum 2 buffers
  1486.  *****************************************************************/
  1487.     void 
  1488. SumFloats(float *in1,  float *in2, float *out, int length)
  1489. /*  in1, in2    ptrs to input buffers
  1490.     out        ptr to output buffer
  1491.     length    length of buffers
  1492. */
  1493. {
  1494. register int    i;
  1495.  
  1496. for (i = 0; i < length; i++) 
  1497.     out[i] = in1[i] + in2[i];
  1498. }    /* ------------------ end SumFloats() --------------- */
  1499.  
  1500. /**************************************************************** 
  1501.  * CopyFloats    copy FLOATs from in buffer to out buffer
  1502.  *****************************************************************/
  1503.     void 
  1504. CopyFloats(float *in, float *out, int length)
  1505. /* in            ptr to input data
  1506.    out            ptr to output data
  1507.    length        length of data */
  1508. {
  1509. int     i;            
  1510.  
  1511. for (i = 0; i < length; i++) 
  1512.     out[i] = in[i];
  1513. }    /* ------------------ end CopyFloats() --------------- */
  1514.  
  1515. /* **********************************************************************
  1516.  * DelayFloats:   delay signal by N samples
  1517.  * ********************************************************************** */
  1518.     void
  1519. DelayFloats(float *in, float *out, long length, 
  1520.         float *delayBuffer, int delayLength, int *delayBufferPtr,
  1521.             char useCircularBuffer)
  1522. /*  in            ptr to input buffer
  1523.     out            ptr to output buffer
  1524.     length        length of input buffer 
  1525.     delayBuffer        ptr to delay buffer
  1526.     delayLength        length of delay buffer
  1527.     delayBufferPtr    index of next read/write position in delay buffer.
  1528.             Need this if delayLength > length
  1529.     useCircularBuffer    if TRUE, use circular buffer to preserve time in
  1530.             sample buffer.  (Usually don't need this.)
  1531. */
  1532. {
  1533. long        i;
  1534. long        delayIndex, delayEndIndex;
  1535. int        firstHalfLength, secondHalfLength;
  1536.  
  1537. /* do delay process w/circular buffer to preserve time in delay buffer */
  1538. if ((delayLength > length)||(useCircularBuffer == TRUE))
  1539.     {
  1540.     delayEndIndex = delayLength-1;
  1541.  
  1542. /* delay line ptr not near delay line end */
  1543.     if (*delayBufferPtr < delayEndIndex-length)
  1544.     {
  1545.     /* write delay buffer into output buffer */
  1546.     CopyFloats(&delayBuffer[*delayBufferPtr], out, length); 
  1547.  
  1548.     /* write input into delay buffer into output buffer */
  1549.     CopyFloats(in, &delayBuffer[*delayBufferPtr], length); 
  1550.  
  1551.     /* advance delay index (bounding check not needed) */
  1552.     *delayBufferPtr += length;
  1553.     }
  1554.  
  1555. /* delay line ptr near delay line end */
  1556.     else
  1557.     {
  1558.     /* write delay buffer into output buffer */
  1559.     firstHalfLength = delayEndIndex - *delayBufferPtr;
  1560.     secondHalfLength = length-firstHalfLength;
  1561.     CopyFloats(&delayBuffer[*delayBufferPtr], out, firstHalfLength);
  1562.     CopyFloats(delayBuffer, &out[firstHalfLength], secondHalfLength);
  1563.  
  1564.     /* write input into delay buffer */
  1565.     CopyFloats(in, &delayBuffer[*delayBufferPtr], firstHalfLength);
  1566.     CopyFloats(&in[firstHalfLength], delayBuffer, secondHalfLength);
  1567.  
  1568.     /* advance and bound delay index */
  1569.     *delayBufferPtr += length;
  1570.     if (*delayBufferPtr > delayEndIndex)
  1571.         *delayBufferPtr -= delayEndIndex;
  1572.     }
  1573.     }
  1574.  
  1575. /* do delay process w/gross copies.  No need to adjust delay buffer ptr.  This 
  1576.     type of delay process doesn't preserve time of samples in delay buffer */
  1577. else
  1578.     {
  1579. /* write delay buffer into output buffer; write input into output buffer */
  1580.     CopyFloats(delayBuffer, out, delayLength); 
  1581.     CopyFloats(in, &out[delayLength], length-delayLength); 
  1582.  
  1583. /* write input tail into delay buffer */
  1584.     CopyFloats(&in[length-delayLength], delayBuffer, delayLength); 
  1585.     }
  1586. }    /* ----------------------- end DelayFloats() --------------------- */
  1587.  
  1588. /****************************************************************
  1589.  * InterleaveNFloatsToShorts interleave N buffers
  1590.  *
  1591.  *            Option to saturate to 16-bits.
  1592.  *                
  1593.  *            Input buffers are assumed to be of equal length
  1594.  *            and data type.
  1595.  *            Output buffer is N*length of input buffer.  First 
  1596.  *            word of input buffer[0] will be first word of 
  1597.  *            output buffer
  1598.  *****************************************************************/
  1599.     void 
  1600. InterleaveNFloatsToShorts(float *inputs[], short *output, int inputLength, 
  1601.                 int interleaveFactor, char saturate)
  1602. /* inputs        ptr to array of input buffers
  1603.    output        ptr to output buffer
  1604.    inputLength        length of input buffer 
  1605.     interleaveFactor    # input buffers 
  1606.    saturate        if TRUE, do saturation 
  1607. */
  1608. {
  1609. register int        i, j;
  1610. register short        *out;
  1611. register float        fValue;
  1612. register float        fCeiling, fFloor;
  1613. register short        iCeiling, iFloor;
  1614.  
  1615. /* rounding = truncation + 1/2 bit DC offset -->> just truncate */
  1616. /* convert buffer, no saturation */
  1617. out = output;
  1618. if (saturate == FALSE)
  1619.     {
  1620.     for (i = 0; i < inputLength; i++) 
  1621.     {
  1622.     for (j = 0; j < interleaveFactor; j++) 
  1623.         out[j] = (short) inputs[j][i];
  1624.     out += interleaveFactor;
  1625.     }
  1626.     }
  1627. /* convert buffer w/saturation to integer range [iFloor..iCeiling] */
  1628. else 
  1629.     {
  1630.     fCeiling =  32767;
  1631.     fFloor   = -32768;
  1632.     iCeiling =  32767;    /* 2^15 - 1 */
  1633.     iFloor   = -32768;    /* -2^15 */
  1634.  
  1635.     for (i = 0; i < inputLength; i++) 
  1636.     {
  1637.     for (j = 0; j < interleaveFactor; j++) 
  1638.         {
  1639.         fValue = inputs[j][i];
  1640.         if        (fValue > fCeiling)
  1641.         out[j] = iCeiling;
  1642.         else if (fValue < fFloor)
  1643.         out[j] = iFloor;
  1644.         else 
  1645.         out[j] = (short) fValue;
  1646.         }
  1647.     out += interleaveFactor;
  1648.     }
  1649.     }
  1650. }    /* ------------------ end InterleaveNFloatsToShorts() --------------- */
  1651.  
  1652. /* **********************************************************************
  1653.  * SecondsToSamples    return time in units of samples
  1654.  * ********************************************************************** */
  1655.     int
  1656. SecondsToSamples(double time, double samplingRate)
  1657. /* time        time (in seconds)  
  1658.   samplingRate    sampling rate (in samples per second) 
  1659. */
  1660. {
  1661. /* compute time in samples
  1662.  *    Seconds * (samples/Second) = samples
  1663.  */
  1664. if ((time >= 0)&&(samplingRate > 0))
  1665.     {
  1666.     return (int) (time*samplingRate);
  1667.     }
  1668. /* report bogus parameters */
  1669. if (time < 0)
  1670.     fprintf(stderr, "SecondsToSamples() bogus time (in seconds) %g\n", time);
  1671. if (samplingRate <= 0)
  1672.     fprintf(stderr, "SecondsToSamples() bogus sampling rate %g\n", samplingRate);
  1673.  
  1674. return (-1);
  1675. }    /* ------------------ end SecondsToSamples() --------------- */
  1676.  
  1677. /* 
  1678.  * LowerPrime     given integer 'number', return prime <= 'number'.
  1679.  *            If 'number' is prime, next prime lesser in value 
  1680.  *            will be returned.  If no lesser prime exists,
  1681.  *            return value will be 2, the least valued prime or
  1682.  *            
  1683.  * ---------------------------------------------------------------------- */
  1684.     int 
  1685. LowerPrime(int number, int lowerBound)
  1686. /* number        number (> 1) to be checked for prime.
  1687.    lowerBound        lowest number (> 1) to check for. (Need not be prime) 
  1688. */
  1689. {
  1690. /* check if input is < 2 */
  1691. if (number < 2)
  1692.     return (2);
  1693.  
  1694. /* check if lowerBound is < 2 */
  1695. if (lowerBound < 2)
  1696.     lowerBound = 2;
  1697.  
  1698. /* check that input is greater than lowerBound */
  1699. if (number < lowerBound)
  1700.     return (lowerBound);
  1701.  
  1702. /* spin around, looking for lesser prime */
  1703. while ((IsPrime(number) == FALSE)&&(number >= lowerBound))
  1704.     number--;
  1705.  
  1706. /* decide whether search reached lowerBound */
  1707. if (number == lowerBound)
  1708.     return (lowerBound);
  1709.  
  1710. return (number);
  1711. }    /* ----------------------- end LowerPrime() --------------------- */
  1712.  
  1713. /* 
  1714.  * UpperPrime    given integer 'number', return prime >= 'number'
  1715.  * ---------------------------------------------------------------------- */
  1716.     int 
  1717. UpperPrime(int number, int upperBound)
  1718. /* number        number (> 1) to be checked for prime.
  1719.     upperBound        highest number (> 1) to check for. (Need not be prime.) */
  1720. {
  1721. /* correct input < 2 */
  1722. if (number < 2)
  1723.     number = 2;
  1724.  
  1725. /* check if lowerBound is < 2 */
  1726. if (upperBound < 2)
  1727.     return (2);
  1728.  
  1729. /* check that input is less than upperBound */
  1730. if (number > upperBound)
  1731.     return (upperBound);
  1732.  
  1733. /* spin around, looking for greater prime */
  1734. while ((IsPrime(number) == FALSE)&&(number <= upperBound))
  1735.     number++;
  1736.  
  1737. /* decide whether search reached upperBound */
  1738. if (number == upperBound)
  1739.     return (upperBound);
  1740. else
  1741.     return (number);
  1742. }    /* ----------------------- end UpperPrime() --------------------- */
  1743.  
  1744. /* 
  1745.  * NearestPrime  given integer 'number', return nearest prime.
  1746.  *            if surrounding primes are equidistant, return greater
  1747.  *            prime.  If numer is prime, it is returned.
  1748.  * ---------------------------------------------------------------------- */
  1749.     int 
  1750. NearestPrime(int number, int lowerBound, int upperBound)
  1751. /* number        number (> 1) to be checked for prime.
  1752.     lowerBound        lowest number (> 1) to check for. (Need not be prime) 
  1753.     upperBound        highest number (> 1) to check for. (Need not be prime)
  1754. */
  1755. {
  1756. int    lowerPrime, upperPrime;
  1757.  
  1758. /* check to see if input is prime */
  1759. if (IsPrime(number) == TRUE)
  1760.     return (number);
  1761.  
  1762. /* find next prime greater in value */
  1763. upperPrime = UpperPrime(number, upperBound);
  1764.  
  1765. /* find next prime lesser in value */
  1766. lowerPrime = LowerPrime(number, lowerBound);
  1767.  
  1768. /*fprintf(stderr, "NearestPrime() loDiff=%d hiDiff=%d\n", 
  1769.         number - lowerPrime, upperPrime - number);*/
  1770.  
  1771. /* determine nearest prime */
  1772. if ((upperPrime - number) <= (number - lowerPrime))
  1773.     return (upperPrime);
  1774.  
  1775. return (lowerPrime);
  1776. }    /* ----------------------- end NearestPrime() --------------------- */
  1777.  
  1778. /* 
  1779.  * IsPrime     given integer 'number', detemine if prime or not and
  1780.  *        return flag (FALSE, TRUE)
  1781.  * ---------------------------------------------------------------------- */
  1782.     int 
  1783. IsPrime(int number)
  1784. /* number        value to test for primeness */
  1785. {
  1786. int i;
  1787. int    isPrime;
  1788. int    truncatedSquareRoot;
  1789.  
  1790. /* if number is even, not prime */
  1791. /* all numbers less than 2 can not be prime */
  1792. if ((((number&0x1) == 0)&&(number != 2))||(number < 2))
  1793.     return(FALSE);
  1794.  
  1795. /* A Composite Number a Will Always Posses A Prime Divisor p Satisfying
  1796.  *    p <= sqrt(a)
  1797.  
  1798.  * Fundamental Theorem of Arithmetic
  1799.     (from Elementary Number Theory, 2nd Ed. by David M. Burton, p.48)
  1800.     Every positive integer n > 1 can be expressed as a product of primes;  
  1801.     this representation is unique, apart from the order in which the factors 
  1802.     occur
  1803.  
  1804.  * from Elementary Number Theory, 2nd Ed. by David M. Burton, p.52
  1805.  
  1806.  *     "There is a property of composite numbers which allows us to reduce
  1807.  * materially the necessary computations - but still the above process
  1808.  * remains cumbersome.  If an integer a > 1 is composite, then it may be
  1809.  * written as a = bc, where 1 < b < a and 1 < c < a.  Assuming that b <= c,
  1810.  * we get b^2 <= bc = a and so b <= sqrt(a).  Since b > 1, Fundamental Theorem
  1811.  * of Arithmetic ensures that b has at least one prime factor p.  Then p <= b
  1812.  * <= sqrt(a);  furthermore, because p|b and b|a, it follows that p|a.  The
  1813.  * point is simply this  a composite number a will always possess a prime
  1814.  * divisor p satisfying p <= sqrt(a)."
  1815.  */
  1816. truncatedSquareRoot = (int) sqrt((double) number);
  1817.  
  1818. /* this can definitely be sped up, w/look up table of primes */
  1819. /* if number is not evenly divisible by any number <= truncatedSquareRoot,
  1820.     number is prime */
  1821. for (i = 2; i <= truncatedSquareRoot; i++)
  1822.     {
  1823.     /* if modulus result is zero, number has been cleanly divided */
  1824.     if (number%i == 0)
  1825.     break;
  1826.     }
  1827. /* figure out what broke loop. If result from modulus operation
  1828.  *     is zero, then number is divisible by smaller number. 
  1829.  *    Thus, number is not prime 
  1830.  */
  1831.     if (i > truncatedSquareRoot)
  1832.     return(TRUE);
  1833.  
  1834. return(FALSE);
  1835. }    /* ----------------------- end IsPrime() --------------------- */
  1836.